#!/usr/bin/python
# -*- coding: utf-8 -*-

# Copyright (c) 2011
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
# Author: Jesus Carrero <j.o.carrero@gmail.com>
#

from scipy import c_, r_, diag, tile, dot
from scipy import zeros, linalg
from numpy import sum as npsum

def twod_inner(tria, quad, func):
    """
    Synopsis:
    --------
      computes the inner product of function f
      against the lagrange function which is one
      on the first node and zero anywhere else.

   Inputt:
   -----
    tria .- 3x2 array.
    quad .- quadrature formula
    f    .- function
  """

    weights = quad.get_weights()
    qpoints = quad.get_qpoints()
    p1_hat = 1.0 - npsum(qpoints, axis=1)

    amat = r_[tria[1, :] - tria[0, :], tria[2, :] - tria[0, :]]
    amat = amat.reshape((2, 2))
    eval_nodes = tile(tria[0, :], (qpoints[:, 0].size, 1))
    eval_nodes += dot(qpoints, amat.T)
    funt_at_quads = func(eval_nodes)

    return linalg.det(amat) * dot(p1_hat, funt_at_quads * weights) / 2.


def oned_inner_for_pide(basis, quad, grid, func, curr_time, next_time):
    """ """

    x_coors = grid.get_grid()
    delta_x = grid.getDx()
    nsegs = grid.getNumSegs()

    qpoints = quad.get_qpoints()
    poly_at_quads = basis.values(qpoints)

    weights = quad.get_weights() / 2.
    poly_at_quads = poly_at_quads * weights.T
    qpoints = (1.0 + qpoints) / 2.

    dof = basis.dof()
    rhs = zeros(((dof - 1) * nsegs + 1, 1))
    for (i, delta_xi) in enumerate(delta_x):
        func_at_quads = func(qpoints * delta_xi + x_coors[i], curr_time, next_time)
        rhs[(dof - 1) * i:(dof - 1) * i + dof] += delta_xi * dot(poly_at_quads,
                func_at_quads)
    return rhs



def oned_inner(basis, quad, grid, func, time=0):
    """ compute
        <\phi, f> with respect to the L2 norm.
        \phi are elements of the basis for the finite
        element space
  """

    x_coors = grid.get_grid()
    delta_x = grid.getDx()
    nsegs = grid.getNumSegs()

    qpoints = quad.get_qpoints()
    poly_at_quads = basis.values(qpoints)

    weights = quad.get_weights() / 2.
    poly_at_quads = poly_at_quads * weights.T
    qpoints = (1.0 + qpoints) / 2.

    dof = basis.dof()
    rhs = zeros(((dof - 1) * nsegs + 1, 1))
    for (i, delta_xi) in enumerate(delta_x):
        func_at_quads = func(qpoints * delta_xi + x_coors[i], time)
        rhs[(dof - 1) * i:(dof - 1) * i + dof] += delta_xi * dot(poly_at_quads,
                func_at_quads)

    return rhs


def ProjOnShapeFunctions(basis, quad, grid, func, time=0):
    """ assemble mass matrix """

    delta_x = grid.getDx()
    nsegs = grid.getNumSegs()

    mass_matrix = basis.get_mass_mat()
    dof = basis.dof()
    diag0 = zeros(((dof - 1) * nsegs + 1, 1))

    sdiag = diag(mass_matrix)
    sdiag = sdiag.reshape((dof, 1))
    for (i, delta_xi) in enumerate(delta_x):
        diag0[(dof - 1) * i:(dof - 1) * i + dof] += sdiag * delta_xi

    diag0 = diag0.flatten()
    if 2 == dof:
        diag1 = mass_matrix[0, 1] * delta_x.flatten()
        mass = diag(diag0, 0) + diag(diag1, 1) + diag(diag1, -1)

    if 3 == dof:
        diag1 = (diag(mass_matrix, 1) * c_[delta_x, delta_x]).flatten()
        diag2 = mass_matrix[0, 2] * c_[delta_x, zeros(delta_x.shape)].flatten()
        diag2 = diag2[0:-1]

        mass = diag(diag0) + diag(diag1, 1) + diag(diag1, -1) \
            + diag(diag2, 2) + diag(diag2, -2)

    rhs = oned_inner(basis, quad, grid, func, time)
    return linalg.solve(mass, rhs)


